Práctica 3
Ayudando a Mathematica a resolver ecuaciones diferenciales ordinarias.

Introducción

En esta práctica te voy a proponer que hagas algunos ejercicios de resolución de EDs sencillas con ayuda de Mathematica. Como hemos comentado en el curso, son pocas las EDs que pueden resolverse, de forma más o menos explícita, calculando primitivas. Algunos tipos de EDs para los que esto es posible son las EDs de variables separadas, las EDs exactas y las que pueden reducirse a ellas como, por ejemplo, la EDO lineal de orden 1 o las EDs homogéneas. Además de estos tipos básicos, hay una gran variedad de EDs cuya solución puede calcularse con técnicas específicas de cada caso. Por ejemplo, las EDs de Bernouilli, de Claireaut, de Ricatti y otras muchas. Saber usar en cada ocasión el recurso adecuado que nos lleve a la solución buscada solamente se aprende con la práctica y, como no se trata de aplicar una teoría sistemática sino una colección de técnicas muy variadas, es algo que se olvida con facilidad. Hay una excepción que es el caso de las EDs lineales y de los sistemas de EDs lineales con coeficientes constantes para los que hay una teoría elaborada que permite, al menos en teoría, obtener de forma explícita sus soluciones.
El propósito de esta práctica es el siguiente. Se trata de que resuelvas algunas EDs con la ayuda de Mathematica. Te voy a proponer algunos ejercicios de resolución de EDs en los que tú deberás reconocer el tipo de ED de que se trata en cada caso y pedirle a Mathematica que realice las operaciones oportunas para obtener la solución. Claro está que Mathematica dispone de un comando que resuelve muchas EDO de los tipos usuales, este comando es DSolve[ ] y está bien documentado en la ayuda del programa. Pero no debes usar este comando para realizar los ejercicios; en todo caso puedes usarlo para comprobar la solución que has obtenido. Es decir, se trata de que tú indiques a Mathematica el procedimiento a seguir y pedirle al programa que haga el "trabajo sucio". De esta forma puedes centrarte en el procedimiento y comprobar su eficacia y puedes descuidarte de cálculos que suelen ser largos. En definitiva, se trata de que tú pongas el cerebro pensante y Mathematica la fuerza bruta. A continuación desarrollo algunos ejemplos para que te sirvan de orientación de lo que pretendo que hagas en esta práctica.
Advertencia importante. Esta práctica debe hacerse con Mathematica 5.1 o posterior.

Dinámica de poblaciones

La ecuación logística

En los llamados procesos de crecimiento o decrecimento exponencial (o maltusiano) se supone que el crecimiento de una población es proporcional al número de individuos. Si representamos por y(t) el número de individuos de la población en el tiempo t, se supone que y' (t) = α y(t) donde α es una constante llamada tasa relativa de crecimiento. La solución de esta sencilla ED es y(t)=(Subscript[y, 0])(e^(α t)) donde (Subscript[y, 0]) es el valor inicial y(0). Si α > 0 la población crecerá exponencialmente y si α < 0 decrecerá exponencialmente. En condiciones ideales este modelo se ajusta, al menos durante cortos perídos de tiempo, al estudio de algunas poblaciones. Pero es claro que el modelo no puede ser válido fuera de esas condiciones ideales y deben tenerse en cuenta factores limitativos del crecimiento que dependerán de los recursos disponibles, depredadores y otras circunstancias que habrá que considerar en cada caso. Una hipótesis más real consiste en suponer que la tasa de crecimiento relativa no es constante sino que depende del número de individuos. Esto lleva a la ED y' (t) = y(t) h(y). Es razonable suponer que h(y) ≃ α > 0 cuando y es pequeño (lo que implica que y' (t) > 0 y la población aumenta) y que h(y) < 0 cuando y es grande (lo que implica que y' (t) < 0 y la población disminuye). La función más simple que cumple estos requisitos es h(t) = α - λ y. Resulta así la ED y' (t) = y(
              t) (α - λ y(t)) = α y(
              t) (1 - β y(t)) llamada ecuación logística. La constante k = 1/β representa un valor llamado capacidad de contencion de la población. Esta ecuación suele escribirse en la forma:

y' (t) = α y(
          t) (1 -  ((y(t))/k))

Busquemos soluciones constantes, y = c de dicha ED. Para ello debe verficarse que α c (1 - c/k) = 0 lo que proporciona c = 0 y c = k . Las soluciones y = 0,  y = k se llaman soluciones de equilibrio pues para dichos valores no se produce cambio en la población. Observa que y' (t) > 0 cuando 0 < y(t) < k (la población crece) y que y' (t) < 0 cuando y(t) > k (la población decrece). 
Para calcular la solución con valor inicial y(0) = (Subscript[y, 0]) observamos que se trata de una ED de variables separadas pues puede escribirse supuesto que y(t) ≠ 0,  y(t) ≠ k, como

\(TraditionalForm\`dy\/\(y(1 - \ y/k)\) = α\ d\ t\)

Ahora ponemos a trabajar a Mathematica y le pedimos que haga lo que nosotros haríamos para resolver esta ED. A saber, calcular las correspondientes primitivas, introducir una constante y tratar de expresar y como función de t (lo que no siempre podrá hacerse). Calculamos después el valor de la constante c por la condición de que para t = 0 se obtenga un valor inicial dado por y(0) = p. Se sustituye el valor de c obtenido y recuperamos la solución como función de todos los parámetros que aparecen en ella, α,  k,  p y t, pues de esta forma podremos estudiar la dependencia de la solución respecto de dichos parámetros.

sol[t_] = 
      y /. Solve[Integrate[(1/(y (1 -  y/k))), y] == α t + c, 
              y][[1]] // Simplify
Off[Solve::"ifun"]
c /. Solve[sol[0] == p, c][[1]];
y[α_, k_, p_][t_] = 
      sol[t] /. c -> % // Simplify
On[Solve::"ifun"]

Observa que el denominador de esta expresión es siempre positivo porque suponemos que α, β, y p son números positivos, por lo que la solución obtenida está definida en todo \(TraditionalForm\`\(\(ℝ\^+\)\(.\)\)\) Comprobamos que los cálculos fectuados son correctos.

{y[α, k, (Subscript[y, 0])][0] (* 
        valor inicial de la población *) , 
Simplify[
        D[y[α, k, p][t], 
            t] == α y[α, k, p][
              t] (1 - (( y[α, k, p][t])/k))] (* 
        comprobamos que es soluciónde la ED *) }

Comprobamos ahora las soluciones de equilibrio y = 0,  y = k.

\({\(y[α, k, 0]\)[t] (*\ 
        y = 0\ es\ una\ solución\ de\ equilibrio\ *) , 
\(y[α, k, k]\)[t] // Simplify (*\ 
        y = k\ es\ una\ solución\ de\ equilibrio\ *) }\)

Aunque ambas son soluciones de equilibrio porque dan lugar a una población constante en el tiempo, la situación de equilibrio que reflejan es de muy distinta naturaleza. Podemos estudiar el comportamiento de la población para valores muy grandes de t. En el siguiente cálculo es obligado informar a Mathematica de que α y p son positivos.

Limit[y[α, k, p][t], t -> +∞, 
      Assumptions -> α > 0 && p > 0]

Esto nos dice que la población tiende a estabilizarse en torno a la solución de equilibrio y = k. Veámoslo gráficamente.

Plot[y[ .2, 1000, 20][t], {t, 0, 50}] ;

Plot[y[ .2, 1000, 1010][t], {t, 0, 50}] ;

La siguiente gráfica muestra el comportamiento de las soluciones cuando α =  .2,  k = 1000 para distintos valores de p.

Show[
        Table[Plot[y[ .2, 1000, p][t], {t, 0, 30}, 
            PlotRange -> {{0, 30}, {0, 1200}}, 
            DisplayFunction -> Identity], {p, 50, 1900, 50}], 
        DisplayFunction -> $DisplayFunction];

Vemos que las soluciones con un valor inicial p > 0 se alejan de la solución y = 0 por lo que se dice que y = 0 es una solución de equilibrio inestable. Vemos también que las soluciones con un valor inicial p > 0 se aproximan de la solución y = k por lo que se dice que y = k es una solución de equilibrio asintóticamente estable.

Modelo poblacional con un umbral crítico

Si en la ecuación logística (1) suponemos que α < 0 entonces el comportamiento de las soluciones es muy diferente y el modelo se ajusta a una situación en la que por debajo de un cierto valor límite, que viene dado por la constante k, la población se extingue, y por encima de dicho valor la población crece indefinidamente. El valor de k es el umbral de supervivencia de la especie.

y[-α, k, p][t] // Simplify

El denominador de esta expresión no se anula nunca si p < k en cuyo caso la población se extingue rápidamente debido a la presencia de la exponencial positiva en el denominador. Podemos verlo gráficamente.

Plot[y[- .2, 1000, 990][t], {t, 0, 
          50}];

Si p > k el denominador se anula en un punto dado por:

λ[α_, k_, p_] = 
      t /. Solve[(E^(t α)) (k - p) + p == 
              0, t][[1]]

En este caso la solución solamente tiene sentido en el intervalo [0, λ[α, k, z]] y tiene una asíntota vertical en el punto λ[α, k, z] lo que indica que la población crecerá sin límite. Podemos verlo gráficamente.

Plot[y[- .2, 1000, 1010][t], {t, 
            0, λ[ .2, 1000, 
                1010] -  .2}] ;

Observa que las soluciones de equilibrio siguen siendo las mismas que antes (no dependen de α) pero ahora y = 0 es asintóticamente estable, mientras que y = k es inestable.

Optimización de recursos perecederos

Supongamos una especie que, en su medio natural, responde al modelo logístico (1), sufre pérdidas sistemáticas por capturas. Podemos considerar que se trata de una población de peces cuyas pérdidas son debidas a la actividad pesquera en la zona que habitan. Naturalmente, los pescadores quieren pescar cuanto más mejor, pero saben que si se exceden ponen en peligro la supervivencia del banco de pesca. Un modelo para esta población suponiendo que se pesca semanalmente un número fijo de peces está dado por la siguiente ED.

y' (t) = α y(
            t) (1 -  ((y(t))/k)) - c

Donde la constante c es el número de peces que se pescan semanalmente. Se trata de una pequeña modificación de la ecuación logística. Observa que la ED (2) sigue siendo una ED de variables separadas porque podemos escribirla en la forma:

\(TraditionalForm\`dy\/\(y\ \((1 - \ β\ y)\) - c/α\) = α\ d\ t\)

Donde hemos puesto β = 1/k. Intentemos proceder igual que antes.

Clear["@"]
sol[t_] = 
      y /. Solve[
              Integrate[(1/(y (1 -  β y) - c/α)), 
                  y] == α t + m, y][[1]] // 
        Simplify
m /. Solve[sol[0] == p, m][[1]];
y[α_, β_, c_, p_][t_] = 
      sol[t] /. m -> % // Simplify

Veremos que la advertencia de Mathematica está justificada. Seguidamente, calculamos las soluciones de equilibrio y simplificamos el resultado obtenido.

y /. 
        Solve[y (1 -  β y) - c/α == 0, 
          y];
soleq[α_, β_, c_] = 
      Simplify[%, α > 0] /. (Sqrt[x_*y_]) -> (Sqrt[x]) (Sqrt[y])

Ahora dichas soluciones dependen de todos los parámetros y solamente tienen sentido para el caso de que sean reales, es decir, cuando α - 4 c β ≥ 0. Observa que la solución que ha calculado Mathematica viene dada por la función tangente. Si reflexionas un poco, caerás en la cuenta de que el resultado debe estar muy simplificado pues la primitiva de la función (1/(y (1 -  β y) - c/α)) depende de que las raíces del denominador sean reales o complejas. Dichas raíces son justamente las soluciones de equilibrio que acabamos de calcular. Cuando son reales, la fracción anterior se descompone en dos fracciones de la forma siguiente.

Apart[(1/(α*β (y - λ)*(y - μ))), 
      y]

Las primitivas correspondientes son logaritmos y, al despejar y en función de t, debemos esperar que se obtengan exponenciales. Aparentemente, en la solución obtenida por Mathematica no hay exponenciales; pero, aunque no se muestren, están ahí escondidas. Recuerda que Mathematica trabaja siempre con números complejos y las funciones usuales están definidas para trabajar por defecto con números complejos. La función Tan[ ] no es una excepción y las funciones trigonométricas están relacionadas con la exponencial compleja. Aquí están las exponenciales que no vemos.

TrigToExp[Tan[z]]

Vamos a tratar de simplificar todo lo que se pueda la solución obtenida para el caso en que α - 4  c b > 0 (raíces reales) y para el caso en que α - 4  c β < 0 (raíces complejas). El caso en que α - 4 c β = 0 (una raíz real doble) hay que tratarlo aparte. Pedimos a Mathematica que simplifique dándole toda la información posible.

Simplify[y[α, β, c, p][t], 
      t > 0 && α > 0 && β > 0 && c > 0 && 
        p > 0 && -α + 4 c β < 0]

Observa que ya no figura la función Tan[ ] sino Tanh[ ], y la tangente hiperbólica se define mediante exponenciales reales. Para que Mathematica simplifique más tenemos que decirle lo que queremos que haga.

% /. 
          Tanh[x_ + 
                y_] -> (Tanh[x] + Tanh[y])/(1 + 
                  Tanh[x] Tanh[y]) /. 
        Tanh[ArcTanh[z_]] -> z;
  
% /. Tanh[x_] -> (Exp[2  x] - 1)/(Exp[2  x] + 1) // 
      Simplify

Ya no puede simplificarse más. Definimos la solución obtenida para el caso de raíces reales.

yreal[α_, β_, c_, p_][
        t_] := (((Sqrt[α]) + (((Sqrt[α - 4 c β]) (-1 + 2 p β - (Sqrt[1 - ((4 c β)/α)]) + (E^(t (Sqrt[α (α - 4 c β)]))) (-1 + 2 p β + (Sqrt[1 - ((4 c β)/α)]))))/(1 - 2 p β + (Sqrt[1 - ((4 c β)/α)]) + (E^(t (Sqrt[α (α - 4 c β)]))) (-1 + 2 p β + (Sqrt[1 - ((4 c β)/α)])))))/(2 (Sqrt[α]) β)) /; (α - 
                4 c β >  0 || 
            c <  (α/(4 β)))

Observa que la solución proporcionada por Mathematica para el caso en que α = 4 c β es y=(1/(2  β))  que es la solución de equilibrio corrrespondiente.

sol[t] /. {c -> (α/(4  β))} // Simplify

Debemos, pues, considerar aparte la ED que corresponde al caso c = (α/(4 β)).

Simplify[α y (1 - β y) - (α/(4 β))]

Se trata, pues, de la ED \(TraditionalForm\`y'\  = \ \(\(-α\)\(\ \)\(\((1 - 2  y\ β)\)\^2/4\) \(β\)\(\ \)\)\). Calculemos la solución que verifica y(0) = p.

sol[t_] = 
      y /. Solve[
              Integrate[(1/(((1 - 2 y β))^2)), 
                  y] == -(α/(4  β)) t + m, 
              y][[1]] // Simplify
m /. Solve[sol[0] == p, m][[1]];
yreal2[α_, β_, p_][t_] = 
      sol[t] /. m -> % // Simplify
Clear[sol]

Consideremos ahora el caso de raíces complejas. Intentemos simplificar.

Simplify[y[α, β, c, p][t], 
      t > 0 && α > 0 && β > 0 && c > 0 && 
        p > 0 && -α + 4 c β > 0]

Como era de esperar, no hemos logrado nada. Debemos hacerlo nosotros.

% /. 
          Tan[x_ + 
                y_] -> (Tan[x] + Tan[y])/(1 - Tan[x] Tan[y]) /. 
        Tan[ArcTan[u_]] -> u;
% /. 
        Tan[w_] -> -I*(Exp[2  I w] - 1)/(Exp[2 I w] + 
                  1);
% /. Exp[I z_] -> Cos[z] + I Sin[z] // Simplify

Eso es todo lo que se puede conseguir. Definamos la correspondiente solución.

yimag[α_, β_, c_, p_][
        t_] := (((Sqrt[α]) + (((Sqrt[-α + 4 c β]) ((Sqrt[-(α/(α - 4 c β))]) (-1 + 2 p β) Cos[(1/2) t (Sqrt[-α (α - 4 c β)])] - Sin[(1/2) t (Sqrt[-α (α - 4 c β)])]))/(Cos[(1/2) t (Sqrt[-α (α - 4 c β)])] + (Sqrt[-(α/(α - 4 c β))]) (-1 + 2 p β) Sin[(1/2) t (Sqrt[-α (α - 4 c β)])])))/(2 (Sqrt[α]) β)) /; (α - 4 c β < 0 || 
            c > (α/(4 β)))

En lo que sigue consideraremos que α = (8/100),  β = (1/1000) (con lo que (α/(4 b)) = 20) y estudiaremos la evolución de la población para distintos valores de p y de c.

Caso de dos raíces reales (c < (α/(4 β))).

Sabemos que siempre que α > 4  c β o, lo que es igual, c < α/4  β, hay dos soluciones de equilibrio. Vamos a ver  que la mayor, es asintóticamente estable y la menor es inestable. Las soluciones de equilibrio para α = (8/100),  β = (1/1000),  
      c = 15, vienen dadas por:

{α, β} = {(8/100), (1/1000)};
soleq[α, β, 15] // Simplify

La solución correspondiente es.

yreal[α, β, 15, p][
        t]//Simplify

Si p < 250 la solución no está definida para todo t ≥ 0 porque el denominador se anula pero el numerador se anula antes por lo que la población se extingue en un tiempo finito. Si p > 250 la solución está definida para todo t ≥ 0 y se tiene que:

Limit[%, t -> +∞, 
      Assumptions -> p > 250]

Es decir, la población se estabiliza en torno al valor y = 750. Podemos verlo gráficamente.

Plot[yreal[α, β, 15, 245][t], {t, 0, 
          70}];
Plot[yreal[α, β, 15, 255][t], {t, 0, 300}, 
        PlotRange -> All];
Plot[yreal[α, β, 15, 755][t], {t, 0, 150}, 
        PlotRange -> All];

Caso de una raíz real doble (c = (α/(4  β))).

En este caso hay una única solución de equilibrio cuyo valor es 1/2  β.

Clear[α, β]
soleq[α, β, α/(4 β)] // 
      Simplify
{α, β, c} = {(8/100), 
          (1/1000), (α/(4 β))};
soleq[α, β, c] // Simplify

La solución de la ED es:

yreal2[α, β, p][t] // Simplify

En este caso, si p > 500 la solución está definida para todo t ≥ 0 y la población se estabilizará asintóticamente en torno a dicho valor:

Limit[%105, t -> +∞, 
      Assumptions -> p > 500]

Si p < 500 el denominador se anula en un punto y la solución no está definida para todo t ≥ 0, pero el numerador se anula antes y los peces se agotarán en un tiempo finito. Se trata de una solución de equilibrio semiestable.

Show[
        Table[Plot[yreal2[α, β, p][t], {t, 0, 40}, 
            DisplayFunction -> Identity], {p, 250, 750, 50}], 
        DisplayFunction -> $DisplayFunction];

Caso de raíces imaginarias (c > (α/(4  β))).

Veamos lo que pasa con las soluciones correspondientes a raíces imaginarias, esto es, c > α/4  β.

Clear[c]
c /:  c > 20 := True
yimag[α, β, c, p][t] // 
        Simplify;
Clear[c];
solimag[c_, p_][t_] = FullSimplify[%%, c > 20 && p > 0]

En este caso la población se extingue siempre en un tiempo finito.

Plot[solimag[21, 1900][t], {t, 0, 100}];

En resumen, el máximo número de peces que se pueden pescar viene dado por c = α/4  β siempre que la población inicial de peces sea mayor que 1/2  β. Si la población inicial de peces es menor que 1/2  β deberá elegirse c < α/4  β con lo cual la población crecerá hasta el punto de equilibrio:

Clear[α, β, c]
soleq[α, β, c][[1]]

el cual es, evidentemente, mayor que 1/2  β. Cuando la población llegue a ser mayor que 1/2  β podrá establecerse como número de capturas  c = α/4  β.

Estudio cualitativo de la naturaleza de los puntos críticos de una ED autónoma

Las EDs del tipo y' = f(y) donde f es una función de ua variable definida en un cierto intervalo, se llaman EDs autónomas. La ED logística es un ejemplo de ellas. Las soluciones de la ecuación f(y) = 0 se llaman puntos críticos o soluciones de equilibrio de la ED. En casos sencillos es posible deducir su naturaleza sin necesidad de resolver la ED. Por ejemplo, en la ED logística se tiene que f(y) = α y (1 - (y/k)), por lo que los puntos críticos son y = 0 e y = k. Representemos gráficamente la curva y = f(y).

<< Graphics`Arrow`
Plot[ .8 y (1 - y/10), {y, 0, 13}, 
        PlotRange -> {- .3, 2}, 
        AxesLabel -> {"y", "f(y)"}, 
        Epilog -> {Text["k/2", {5,  .1}], 
            Text["k", {10.1,  .1}], Arrow[{ .5,  .1}, {3.5,  .1}], 
            Arrow[{6,  .1}, {9,  .1}], 
            Arrow[{13,  .1}, {11,  .1}]}];

Tenemos que f(y) > 0 para 0 < y < k y f(y) < 0 para y > k. Por tanto y' (t) > 0 para 0 < y < k por lo que las soluciones serán crecientes (como funciones de t) cuando y esté en el intervalo [0, k]. Análogamente, las soluciones serán decrecientes para y > k. Las flechas en la gráfica indican esta situación. Deducimos que las soluciones se aproximarán al valor y = k y se alejarán del valor y = 0. Es decir, la solución y = k es estable, mientras que y = 0 es inestable.
Observa que la situación cambia por completo si suponemos que α < 0.

Plot[- .8 y (1 - y/10), {y, 0, 13}, 
        PlotRange -> { .3, -2}, 
        AxesLabel -> {"y", "f(y)"}, 
        Epilog -> {Text["k/2", {5,  .1}], Text["k", {9.9,  .1}], 
            Arrow[{3.5,  .1}, { .5,  .1}], Arrow[{9,  .1}, {6,  .1}], 
            Arrow[{11,  .1}, {13,  .1}]}];

En este caso la solución y = k es inestable, mientras que y = 0 es estable (estamos considerando solamente soluciones positivas, si se permiten soluciones con valores negativos entonces y = 0 sería un punto de equilibrio inestable). Un estudio análogo puede hacerse para la ED y' = α y (1 - β y) - c.

Un criterio fácil para decidir la naturaleza de un punto crítico es el siguiente. Supongamos que y=(Subscript[y, 0]) es un punto cr´tico de la ED y' = f(y). Es decir f((Subscript[y, 0])) = 0. Entonces, si f' ((Subscript[y, 0])) < 0 se tendrá que f' (y) < 0 en un pequeño intervalo [(Subscript[y, 0]) - δ, (Subscript[y, 0]) + δ ] en el cual fserá decreciente. Por tanto f(y) > 0 para y - δ < y < (Subscript[y, 0]), y f(y) < 0 para (Subscript[y, 0]) < y < (Subscript[y, 0]) + δ; luego y' (t) será creciente a la izquierda de (Subscript[y, 0])  y decreciente a la derecha de (Subscript[y, 0]) por lo que (Subscript[y, 0]) será un punto de equilibrio estable. Análogamente, si f' ((Subscript[y, 0])) > 0 el punto (Subscript[y, 0]) será un punto de equilibrio inestable.

Ejercicio 1

y' (t) = 
      r y(t) (1 -  ((y(t))/K)) - m y(t)

Ejercicio 2

Considera la ED autónoma
					d(y/d)t = (r - s) y - 
              a (y^3)
donde \(TraditionalForm\`\(\(r\)\(,\)\(\ \)\(s\)\(,\)\(\ \)\(a\)\(\ \)\)\)son constantes positivas.

Poligonales de Euler

Consideremos una ED de la forma y' = f (x, y). Podemos representar gráficamente las curvas integrales de esta ED sin más que tener en cuenta que la pendiente en un punto (a, b) de una curva integral que pase por  dicho punto tiene que ser igual a f(a, b). Por tanto, podemos aproximar la gráfica de dicha curva por una poligonal, llamada una poligonal de Euler, cuyos vértices vienen dados por

((Subscript[x, 1]), (Subscript[y, 1])) = (a, b); \(TraditionalForm\`\(\(\((x\_2, y\_2)\) = \((a, b)\) + 
          h\ \((1, f(a, b))\)\/\(\(||\)\((1, f(a, b))\)\(||\)\); \ \((x\_3, 
          y\_3)\) = \((x\_2, y\_2)\) + 
          h \(\(\ \)\((1, f(x\_2, y\_2))\)\)\/\(\(||\)\((1, f(x\_2, y\_2))\)\(||\)\)\)\(,\)\(\(.\)\(\ \)\(.\)\(\ \)\(.\)\)\(,\)\(\ \)\)\)\(TraditionalForm\`\(\(\ \)\(\((x\_\(n + 1\), 
          y\_\(n + 1\))\) = \((x\_n, y\_n)\) + 
          h\ \((1, f(x\_n, y\_n))\)\/\(\(||\)\((1, f(x\_n, y\_n))\)\(||\)\)\)\)\), . . .

Es decir, nos desplazamos de un punto (x, y) al siguiente siguiendo la dirección del vector tangente  calculado en dicho punto y normalizado, (((1, 
                f(x, y)) )/(||(1, f(x, y))||)), yavanzando una distancia igual a |h| ( h puede ser una cantidad positiva o negativa). Naturalmente, es de esperar que cuanto más pequeño sea el número h mejor será la aproximación obtenida. Lógicamente, la gráfica debe representarse en algún intervalo que contenga al punto a, digamos [c, d] con c ≤ a ≤ d. La forma de proceder será calcular la poligonal a partir del punto (a, b) hasta llegar al valor de abscisa x = c y unirla apropiadamente a la poligonal que parte de  (a, b) hasta llegar al valor de abscisa x = d. Para ello dividiremos cada uno de los intervalos [c,a] y [a,d] en n partes de igual longitud, definiremos k = (c - a)/n ≤  0 y h = (d - a)/n ≥  0 e iteraremos el proceso antes descrito para obtener los puntos de la poligonal de Euler.

Ejercicio 3

Ecuaciones diferenciales exactas

Te recuerdo que un campo vectorial F(x, y)=(P(x, y), 
          Q(x, y)) , que suponemos definido en un dominio simplemente conexo del plano, es conservativo cuando ∂(Q/∂)x (x, 
            y) = ∂(P/∂)y (x, y) . En tal caso se verifica que hay una función U(x, y), llamada función potencial, única salvo constantes aditivas, que verifica que ∂(U/∂)x = 
        P,  ∂(U/∂)y = Q . Cuando esto es así, se dice que la ED P(x, y) dx +  Q(x, 
              y) dy = 0 es exacta. La solución general de dicha ED es la familia de curvas definidas implícitamente por U(x, y) = C donde C es una constante.

Ejercicio 4

Simplify[x (x + y) - (x^2) - x y == 0]   
Simplify[
      x (x + y) - (x^2) == 
        0]  
Simplify[
      x (x + y) - (x^2) - x y === 
        0]  (* ¡salida False! Por eso no se debe usar  ===  
        como control boleano cuando trabajamos simbólicamente *)  (* combinando los dos comandos se puede lograr algo mejor que parece que funciona bien *) 
Simplify[x (x + y) - (x^2) == 0] === True   
Simplify[
        x (x + y) - (x^2) - x y == 0]===True

Ecuaciones diferenciales lineales

Resolución de una ED lineal homogénea con coeficientes constantes

Raíces reales simples

Sea la ED

y''' - 6  y'' + 11  y' - 6  y = 0

En la siguiente celda definimos la ecuación, calculamos las raíces del polinomio característico, comprobamos si todas son distintas, obtenemos la solución general de la ED y la comprobamos. Observa lo fácil que es, cuando no hay raíces múltiples como ocurre en este caso, definir la solución de la ED aprovechando que en Mathematica casi todas las funciones pueden evaluarse en listas.

Clear["@"]
eqh[y_, x_] = y'''[x] - 6  y''[x] + 11  y'[x] - 
          6  y[x];
Print["Las raíces del polinomio caracterítico \
son:"]
raices = 
      z /. Solve[z^3 - 6  z^2 + 11  z - 6 == 0, 
          z]
Print["La solución general de la ED ",  
      eqh[y, x],  "=0 es:"]
solh[x_] = 
      Exp[raices*x] . {(Subscript[C, 1]), (Subscript[C, 2]), (Subscript[C, 3])}
Print["Comprobación."]
Simplify[eqh[solh, x] == 0]

Raíces reales e imaginarias simples

Vamos a ver una forma de obtener la solución para el caso en que haya raíces imaginarias y reales simples. Consideremos la ED

y''' - 3  y'' + 4  y' - 2  y = 0

Clear["@"]
eqh[y_, x_] = y'''[x] - 3  y''[x] + 4  y'[x] - 
          2  y[x];
Print["Las raíces del polinomio caracterítico \
son:"];
raices = 
      z /. Solve[-2 + 4 z - 3 (z^2) + (z^3) == 0, 
          z]
Print["Calculamos las exponenciales."];
Exp[raices*x]
Print["Formamos las soluciones asociadas a cada \
raíz."];
ComplexExpand[Re[%%]] // Union
% /. {Cos[x_] -> {Cos[x], Sin[x]}, 
          Exp[u_]*Cos[x_] -> {Exp[u]*Cos[x], Exp[u]*Sin[x]}} // 
      Flatten
Print["La solución general de la ED ",  
      eqh[y, x],  "=0 es:"]
solh[x_] = %% . {(Subscript[C, 1]), (Subscript[C, 2]), (Subscript[C, 3])}
Print["Comprobación."]
Simplify[eqh[solh, x] == 0]

Raíces múltiples reales

Consideramos ahora el caso en que hay raíces múltiples. Sea la ED

y''''' + y'''' - 4  y''' - 4  y'' + 11  y' - 5  y = 
      0

Clear["@"]
eqh[y_, x_] = y'''''[x] + y''''[x] - 4  y'''[x] - 
          4  y''[x] + 11  y'[x] - 
          5  y[x];
Print["Las raíces del polinomio caracterítico \
son:"];
raices = 
      x /. Solve[-5 + 11 x - 4 (x^2) - 4 (x^3) + (x^4) + (x^5) == 0, 
          x]
Print["Agrupamos las raíces"];
Split[raices]
Print["Calculamos las exponenciales"];
Exp[%%*x]
Print["Tenemos en cuenta las soluciones correspondientes a ra\
íces múltiples"];
%% /. z_?VectorQ :> z*Table[x^k, {k, 0, Length[z] - 1}] // 
      Flatten
Print["Formamos las soluciones asociadas a cada \
raíz."];
ComplexExpand[Re[%%]] // Union
Print["Completamos la solución"];
%% /. {Cos[x_] -> {Cos[x], Sin[x]}, 
          Exp[u_]*Cos[x_] -> {Exp[u]*Cos[x], Exp[u]*Sin[x]}} // 
      Flatten
Print["La solución general de la ED ",  
        eqh[y, x],  "=0 es:"];
solh[x_] = {(Subscript[C, 1]), (Subscript[C, 2]), (Subscript[C, 3]), (Subscript[C, 4]), (Subscript[C, 5])} . %%
Print["Comprobación."];
Simplify[eqh[solh, x] == 0]

Raíces reales e imaginarias múltiples y simples

Consideremos la ED lineal de orden 11:

32  y[x]-112  y'[x]+184 y''[x]-180 Derivative[3][y][x]+80 Derivative[4][y][x]+68 Derivative[5][y][x]-174 Derivative[6][y][x]+179 Derivative[7][y][x]-112 Derivative[8][y][x]+44 Derivative[9][y][x]-10 Derivative[10][y][x]+Derivative[11][y][x]=0

Veamos cómo obtener el polinomio característico a partir de la ED.

Clear["@"]
eq[y_, x_]=32  y[x]-112  y'[x]+184 y''[x]-180 Derivative[3][y][x]+80 Derivative[4][y][x]+68 Derivative[5][y][x]-174 Derivative[6][y][x]+179 Derivative[7][y][x]-112 Derivative[8][y][x]+44 Derivative[9][y][x]-10 Derivative[10][y][x]+Derivative[11][y][x];
policar[
        x_] = eq[y, x] /. Derivative[k_][y][x] -> (x^k) /. 
        y[x] -> 1

Podemos resolver esta ED homogénea aplicando lo que ya sabemos.

Clear["@"]
eqh[y_, x_]=32  y[x]-112  y'[x]+184 y''[x]-180 Derivative[3][y][x]+80 Derivative[4][y][x]+68 Derivative[5][y][x]-174 Derivative[6][y][x]+179 Derivative[7][y][x]-112 Derivative[8][y][x]+44 Derivative[9][y][x]-10 Derivative[10][y][x]+Derivative[11][y][x];
  
policar[
          x_] = eqh[y, x] /. Derivative[k_][y][x] -> (x^k) /. 
          y[x] -> 
            1;
Print["Las raíces del \
polinomio caracterítico son:"];
raices = x /. Solve[policar[x] == 0, 
          x]
Print["Agrupamos las \
raíces"];
Split[
      raices]
Print["Calculamos las \
exponenciales"];
Exp[%%*
        x]
Print["Tenemos en cuenta las soluciones \
correspondientes a raíces múltiples"];
%% /. 
        z_?VectorQ :> z*Table[x^k, {k, 0, Length[z] - 1}] // 
      Flatten
Print["Calculamos la parte real y eliminamos \
elementos repetidos y ordenamos"];
ComplexExpand[Re[%%]] // 
      Union
Print["Completamos la solución"];
%% /. {Cos[x_] -> {Cos[x], Sin[x]}, 
          Exp[u_]*Cos[x_] -> {Exp[u]*Cos[x], Exp[u]*Sin[x]}} // 
      Flatten
Print["La solución general de la ED ",  
        eqh[y, x],  "=0 es:"];
solh[x_] = %% . 
        Table[
          (Subscript[C, k]), {k, 1, 
            Length[raices]}]
Print["Comprobación."];
Simplify[
      eqh[solh, x] == 0]

Ejercicio 5

Cálculo de soluciones particulares

Método de los coeficientes indeterminados

Ejemplo 1

Calculemos una solución particular de la ED

y''' - 5  y'' + 8  y' - 4  y = 
      cos (4  x)

El término independiente es de la forma \(TraditionalForm\`cos(
          4  x) = \(P(
            x)\) \(e\^\(α\ x\)\) \(cos(β\ x)\)\)  por lo que α = 0,  β = 4,  
      P(x) = 1. Tenemos que α + i β = 
        4  i no es raíz del polinomio característico. Buscaremos una solución particular de la forma y(x)=z (e^(4  i x)) donde z es una constante que puede tomar valores complejos que calculamos por la condición de que y(x) sea solución de la ED  \(TraditionalForm\`y''' - 5  y'' + 8  y' - 
          4  y = \(\(P(
              x)\) e\^\(\((α + i\ β)\) x\) = e\^\(4  i\ x\)\)\). Finalmente, la solución que necesitamos será la parte real de la función y(x)

Clear["@"]
eq[y_, x_] := y'''[x] - 5  y''[x] + 8  y'[x] - 
          4  y[x];
y[x_] = z*Exp[4 I x];
eq[y, x] // Simplify

Imponemos que y sea solución de la ED y''' - 5  y'' + 8  y' - 
          4  y = (e^(4  i x)).

z = z /. Solve[% == Exp[4 I x], z][[1]]

Una vez obtenido z calculamos la parte real de la solución obtenida y comprobamos el resultado.

h[x_] = ComplexExpand[Re[y[x]]]
eq[h, x] // Simplify

Ejemplo 2

Busquemos ahora una solución particular de la ED

y'''' - 5  y'' + 8  y' - 4  y = 
      x sen(x)

El término independiente es de la forma \(TraditionalForm\`sen(
          x) = \(P(
            x)\) \(e\^\(α\ x\)\) \(sen(β\ x)\)\)  por lo que α = 0,  β = 1,  
      P(x) = 1. Tenemos que \(TraditionalForm\`α + i\ β = i\) no es raíz del polinomio característico. Buscaremos una solución particular de la forma y(x)=(z + w x) (e^(i x)) donde z,  w son constantes que pueden tomar valores complejos que calculamos por la condición de que y(x) sea solución de la ED  \(TraditionalForm\`y''' - 5  y'' + 8  y' - 
          4  y = \(\(P(
              x)\) e\^\(\((α + i\ β)\) x\) = x\ e\^\(i\ x\)\)\). Finalmente, la solución que necesitamos será la parte imaginaria de la función y(x)

Clear["@"]
eq[y_, x_] := y'''[x] - 5  y''[x] + 8  y'[x] - 
          4  y[x];
y[x_] = (z + w x) Exp[I x];
eq[y, x] // Simplify

Observa que siempre que se sustituye en una EDL homogénea una función de la forma \(TraditionalForm\`\(Q(
          x)\) e\^\(\((α + i\ β)\) x\)\) podemos sacar factor común la función exponencial. Es conveniente eliminar esa función para identificat coeficientes.

CoefficientList[Exp[-I x] %, x]

Ahora identificamos estos coeficientes con los de P(x) = x.

{z, 
        w} = {z, 
          w} /. Solve[% == {0, 1}, {z, w}][[1]]

Una vez obtenidos z,  w calculamos la parte imaginaria de la solución obtenida y comprobamos el resultado.

h[x_] = ComplexExpand[Im[y[x]]]
eq[h, x] // Simplify

Ejemplo 3

Resolvemos la ED

y''''' + y'''' - 4  y''' - 4  y'' + 11  y' - 
        5  y = (1 - x) (e^(-2 x)) cos(
          x)

Calculamos las raíces del polinomio característico.

x /. Solve[-5 + 11 x - 4 (x^2) - 4 (x^3) + (x^4) + (x^5) == 0, 
        x]

El término independiente es de la forma \(TraditionalForm\`x\ \(e\^\(\(-2\) x\)\) \(cos(
            x)\) = \(P(
            x)\) \(e\^\(α\ x\)\) \(cos(β\ x)\)\)  por lo que α = -2,  β = 1,  
      P(x) = 1 - x. Tenemos que \(TraditionalForm\`α + i\ β = \(-2\) + i\)  es raíz simple del polinomio característico. Buscaremos una solución particular de la forma y(x)=x(z +  w x) \(e\^\(\((\(-2\)\  + \ i)\)\ x\)\) donde z,  w son constantes que pueden tomar valores complejos que calculamos por la condición de que y(x) sea solución de la ED   y'''' + 5  y'' + 4  y=P(x) \(e\^\(\((\(-2\)\  + \ i)\)\ x\)\). Finalmente, la solución que necesitamos será la parte real de la función y(x) .

Clear["@"]
eq[y_, x_] = y'''''[x] + y''''[x] - 4  y'''[x] - 
          4  y''[x] + 11  y'[x] - 
          5  y[x];
y[x_] = 
        x (z + w x) 
          Exp[(-2 + I) x];
eq[y, x] // Simplify

CoefficientList[Exp[(2 - I) x]*%, x]

Identificamos coeficientes.

{z, 
        w} = {z, 
          w} /. Solve[% == {1, -1}, {z, w}][[1]]

Una vez obtenidos z,  w calculamos la parte real de la solución obtenida y comprobamos el resultado.

h[x_] = ComplexExpand[Re[y[x]]]
eq[h, x] // Simplify

Ejercicio 6

Método de variación de constantes

El método de los coeficientes indeterminados solamente puede aplicarse cuando el término independiente es de la forma \(TraditionalForm\`\(P(
          x)\) \(e\^\(\(α\)\(\ \)\(x\)\(\ \)\)\) \(cos(β\ x)\)\) o \(TraditionalForm\`\(P(
          x)\) \(e\^\(\(α\)\(\ \)\(x\)\(\ \)\)\) \(sen(β\ x)\)\) o una suma de ellos. El método que vamos a considerar ahora es más general porque puede aplicarse siempre. Veamos un ejemplo.

Consideremos la ED

y''' - 3  y'' + 4  y' - 
        2  y = (e^x) tg(x)

Lo primero que debemos hacer es calcular la solución general de la EDL homogénea correspondiente. Dicha solución es

\(TraditionalForm\`\(C\_1\) e\^x\  + \(C\_2\) e\^x\ Cos[x]\  + \(\(C\_3\) \(e\^x\)\(\ \)\(Sin[
          x]\)\(\ \)\)\)

En el método de variación de constantes se busca una solución particular de la forma

(Subscript[C, 1])(x) (e^x)  + (Subscript[C, 2])(
          x) (e^x) Cos[x]  + (Subscript[C, 3])(
          x) (e^x) Sin[x]

donde (Subscript[C, 1])(x), (Subscript[C, 2])(x) y (Subscript[C, 3])(x) son funciones que hay que calcular. Para ello se procede siempre de la misma forma.

Definimos una función auxiliar como sigue.

aux[
          x_] = (Subscript[a, 1]) Exp[x] + (Subscript[a, 2]) Exp[x] Cos[x] + (Subscript[a, 3]) Exp[x] 
            Sin[x];

Estas funciones (Subscript[a, 1])(x),  (Subscript[a, 2])(x),  (Subscript[a, 3])(x) representan las derivadas de las funciones  (Subscript[C, 1])(x), (Subscript[C, 2])(x) y (Subscript[C, 3])(x) que queremos calcular. Como no nos interesa que Mathematica las derive las hemos representado como si fueran constantes. Ahora resolvemos el sistema de n ecuaciones (tantas como el orden de la ED) lineales en las incógnitas (Subscript[a, 1])(x),  (Subscript[a, 2])(x),  (Subscript[a, 3])(x) que se obtiene igualando la función "aux" y sus derivadas sucesivas hasta la de orden n-2 a cero, y la derivada de orden n-1 se iguala al término independiente de la ED. Este sistema permite calcular las funciones  (Subscript[a, 1])(x),  (Subscript[a, 2])(x),  (Subscript[a, 3])(x).

{(Subscript[a, 1]), (Subscript[a, 2]), 
        (Subscript[a, 3])} /. Solve[{aux[x], aux'[x], aux''[x]} == {0, 0, 
              Exp[x] Tan[x]}, {(Subscript[a, 1]), (Subscript[a, 2]), (Subscript[a, 3])}][[1]]

Ahora calculamos las primitivas de las funciones obtenidas.

Integrate[%, x]

La solución particular buscada es la siguiente.

solpart[x_] = % . {Exp[x], Exp[x] Cos[x], Exp[x] Sin[x]} // 
        Simplify

Lo comprobamos sustituyendo esta función en la ED.

eqd[y_, x_] = y'''[x] - 3  y''[x] + 4  y'[x] - 
          2  y[x];
eqd[solpart, x] // Simplify

Ejercicio 7

y'' - 2  y' + y = (e^x) 
        arctg x